Role of electronic localization in the phosphorescence of iridium sensitizing dyes 
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In this work we present a systematic study of three representative iridium dyes, namely, Ir(ppy)3, 
FIrpic and PQIr, which are commonly used as sensitizers in organic optoelectronic devices. We 
show that electronic correlations play a crucial role in determining the excited-state energies in 
these systems, due to localization of electrons on Ir d orbitals. Electronic localization is captured by 
employing hybrid functionals within time-dependent density-functional theory (TDDFT) and with 
Hubbard-model corrections within the A-SCF approach. The performance of both methods are stud- 
ied comparatively and shown to be in good agreement with experiment. The Hubbard-corrected 
functionals provide further insight into the localization of electrons and on the charge-transfer char- 
acter of excited-states. The gained insight allows us to comment on envisioned functionalization 
strategies to improve the performance of these systems. Complementary discussions on the A-SCF 
method are also presented in order to fill some of the gaps in the literature. 
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I. INTRODUCTION 

Phosphorescent organomctallic dyes have attracted 
considerable interest over the last few decades as 
highly efficient sensitizers in organic light-emitting diodes 
(OLEDs) Q. These complexes are characterized by 
strong spin-orbit coupling due to heavy transition metals 
in their core, such as Ir, Pt, and Os. The strong spin-orbit 
coupling leads to intersystem crossing from the singlet to 
triplet excited states and allows for the emission from 
the otherwise forbidden triplet state (phosphorescence), 
in addition to emission from the singlet state (fluores- 
cence). As a result, phosphorescent emitters can have an 
internal quantum efficiency of 100% [l|. Moreover, the 
possibility to tune the colors of these emitters by modi- 
fications in the ligands surrounding the transition-metal 
center has opened up the possibility to design color dis- 
plays and efficient white light sources from organic ma- 
terials (see Ref. for a review). 

More recently, the strong intersystem-crossing mecha- 
nism that take place in these complexes have been ex- 
ploited to design organic solar cells (OSCs) of improved 
power conversion efficiency [HH]]. In fact, in organic het- 
crojunctions, sensitizing transition-metal complexes are 
utilized as promoters in converting photogenerated sin- 
glet excitons into long-lived triplet excitons, thereby in- 
creasing the probability for an exciton to diffuse to a 
donor-acceptor interface and dissociate into collectible 
charges In a recent study @ it was reported that 

incorporating Ir(ppy)3 (in concentrations as low as few 
mass percents) into electron-donating polymers nearly 
doubled the photovoltaic efficiency of the considered or- 
ganic heterojunctions. 



These remarkable experimental achievements motivate 
further theoretical effort to understand the optical prop- 
erties of transition-metal complexes and ultimately opti- 
mize their sensitizing functions taking advantage of their 
almost unlimited chemical versatility, i.e., modifying the 
transition-metal center and attached chromophores. To 
assist in this endeavor, computational methods of in- 
creasing predictive ability are now available. Among suc- 
cessful electronic-structure techniques, the Green's func- 
tion many-body perturbation theory (GW) with Bethc- 
Salpeter post-treatment has been accurate in predicting 
excited-state energies [11, yet requiring considerable 
effort to achieve full self-consistent convergence. Com- 
putationally less demanding predictions are based upon 
time-dependent density-functional theory (TDDFT) and 
linear- response theory 04^1 ■ Recently, alternative A- 
SCF [l(| techniques that consist in evaluating excited- 
state energies from moderately demanding constrained 
density-functional theory (DFT) calculations have also 
seen a revival of interest }ll| - [l3| . 

In spite of their limited cost, the accuracy of TDDFT 
and A-SCF approximations depends crucially on the 
ability of the underlying DFT functional to properly cap- 
ture orbital localization. In this regard, local and scmilo- 
cal DFT functionals are known to insufficiently localize 
electronic states. Moreover, conventional TDDFT cal- 
culations rely on the adiabatic approximation whereby 
the frequency dependence of the TDDFT kernel is ne- 
glected, representing another important source of error 
in capturing excited states. In practical terms, TDDFT 
approximations do not properly describe charge-transfer 
states that consist of an electron weakly coupled to a 
hole MM- 
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In this study, we examine the performance of TDDFT 
and A-SCF approximations in predicting singlet and 
triplet excitonic energies in three representative Ir com- 
plexes, namely, Ir(ppy)3, FIrpic and PQIr, which emit 
in the green, the blue and the red [22|], respectively. 
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Specifically, wc compare TDDFT with A-SCF that em- 
ploy Hubbard corrections, enabling us to gain needed in- 
sight into the effects of electronic localization and charge- 
transfer for excited states and comment on envisioned 
functionalization strategies. Such insight will be useful 
to study the influence of ligand modification on singlet 
and triplet energies with the ultimate goal of guiding ex- 
periments in designing more efficient OLEDs and OSCs. 

The paper is organized as follows: We provide some 
technical details about the computational approach used 
in this work in Section|ll]and summarize the methods em- 
ployed to compute excited-state energies within TDDFT 
and A-SCF in Section InTl In Section HVl we present the 
main results and provide a detailed discussion on the 
performance of the methods used and some comments 
on strategies to tune excited-state energies. Finally, we 
provide some concluding remarks. We also include an 
Appendix, which provides a brief analysis of the calcula- 
tion of the singlet excited-state energy, within the A-SCF 
approach. 



III. 



THEORETICAL OVERVIEW 



II. COMPUTATIONAL METHODS 



Structural optimizations and A-SCF calculations pre- 
sented in this paper were performed using the plane- 
waves pscudopotcntial implementation of DFT contained 
in the PWSCF code of the Quantum ESPRESSO pack- 
age [23j]. The TDDFT calculations were performed 
using the Gaussian 09 package 24|. For plane- wave 
calculations, the exchange-correlation energy was ap- 
proximated using the generalized-gradient approxima- 
tion (GGA) with the Perdew-Burke-Ernzherof (PBE) 
parametrization (2f|. The Ir, N, C, and H atoms were all 
represented by ultrasoft pseudopotentials [2(| . The elec- 
tronic wavefunctions and charge density were expanded 
up to kinetic energy cutoffs of 30 Ry and 360 Ry, respec- 
tively. In performing TDDFT calculations, the SDD ba- 
sis set of Refs. I27T4331 was used. The semilocal PBE (2f| 
and hybrid functional B3LYP [Mill and M06 [H] were 
used in the TDDFT calculations. 



In our plane wave calculations, we employed the 
Hubbard-model corrected DFT+U [371 - 140} method, in 
order to capture the effects of electronic localization ac- 
curately. An improved version of the DFT+U formalism, 
which includes inter-site interactions (DFT+U+V) [4lj . 
was also used in order to obtain a better description of 
interactions between localized electrons on Ir d states 
and the surrounding organic ligands. The values of the 
Coulomb interaction parameters U and V were computed 
using the linear-response method introduced in Ref. |40| . 
The molecular structures and charge densities presented 
in this paper were generated using XCrysden 42 1. 



In this section, we provide a brief overview of the com- 
putational approaches used in this paper. We refer the 
reader to the original literature for a complete discussion 
of these methods. 



A. TDDFT 

The theorem of Runge and Gross Q states that there 
is a one-to-one correspondence between a time-dependent 
external potential and the electronic density of a many- 
body system. A particular application of the Runge- 
Gross theorem was used to compute the energies of ex- 
cited states of many-body systems, using linear-response 
functions @, Q . Within this approach, the poles of the 
response function provide the excited-state energies of 
the system, which can be determined through the solu- 
tion of the following eigenvalue problem [i| [27l - l33j : 



L M 

M* L* 



X(w) 
Y( W ) 



-I 
I 



Y( W ) 



(1) 



where the solutions uj represent excited-state energies. 
The matrices L and M are given by 

a a;j b r — ^cr r $i j &a b (^a ) ^-i a a\b j r 

M iarT , jbr = I d 3 r [ d 3 r^r(r)€(r) 



xi^(r,r')V>n r ')^ T (r), 



+ 



(Spa(r) Sp T (r') 



(2) 



where ip° are solutions to the Kohn-Sham (KS) equa- 
tions with eigenvalues e° , E xc is the exchange-correlation 
functional, and p a denotes the single-particle density con- 
structed from occupied Kohn-Sham states. In addition, 
we adopt the convention that the indices i, j run over oc- 
cupied states, while a, b run over unoccupied states, and 
the frequency-dependent vector coefficients Xi a a = Y a i a 
are related to the response of the electronic density to a 
time-dependent perturbation through 

5 Pa{r, t) = [du ]T [l,„(,)C(r)<(r)e- lui + (a o »)] 

i. a 

(3) 

Very often, the frequency dependence of the approximate 
exchange-correlation functional E xc is ignored (the adia- 
batic approximation), in obtaining solutions to Eq. ([1]). 
It is well known that the approach briefly summarized 
above usually fails in reproducing experimental excited- 
state energies for Rydberg and charge-transfer-type ex- 
citations [l3 - l2l| . In charge-transfer excitations, the 
occupied ipf and unoccupied ip° orbitals are spatially 
separated. Therefore, the matrix M, representing the 
off-diagonal block in Eq. ([T]) becomes negligible for lo- 
cal or semilocal exchange-correlation functionals (such as 
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PBE) , while L reduces to a diagonal matrix composed of 
energy differences of occupied and unoccupied KS states. 
As a result, the excitation energies are simply given by 
KS energy differences, which underestimate the experi- 
mental values 1J]. This inaccuracy can be corrected to 



some extent with the use of nonlocal, hybrid functionals 
(such as B3LYP and M06). The inaccuracy related to 
Rydberg and charge-transfer type excitations stems from 
the incorrect asymptotic behavior of most approximate 
exchange-correlation functionals, and can be mitigated 
by improving their long-range behavior [l9|, [U l43Tj50| . 
The modifications to the long-range interactions are gen- 
erally parameterized and determined empirically. Alter- 
natively, they can be determined from first-principles in 
a system-dependent fashion Several works also 

considered the A-SCF method (discussed in the next sub- 
section) to be more precise for these types of excitations 
compared to TDDFT [HHlI El- A mixture of the two 
approaches has also been shown to be effective in treating 
charge-transfer excitations [HI, EH . 



B. A-SCF 

A-SCF method is based on the construction of excited- 
state electronic densities, using a non-Aufbau scheme for 
orbital occupations [To| . For a system of N electrons, 
the electronic density is constructed by filling the lowest 
N — 1 orbitals and the (N + l) th orbital at each itera- 
tion of the electronic-structure optimization until self- 
consistency is reached. More precisely, the electronic 
density of an excited-state of spin a is given by 



N-1 

£ 



lVf(r)| 2 + |VWr)| 



(4) 



where the orbitals tpf are self-consistently determined 
from the minimization of the energy functional depend- 
ing on the density in Eq.Q. Examples of such configu- 
rations are schematically represented in Fig.[TJ where the 
ground-state is assumed to be in a closed shell configu- 
ration. The constructed density can be reproduced from 
a single excited Slater determinant of specific space and 
spin symmetry (for some cases). In Fig. [TJ $g refers to 
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FIG. 1: Schematic representation of the single determinant 
configurations used in A-SCF calculations 

the ground-state wavefunction, $t is a triplet excited- 
state, while 'I'mi and $M2 are excited-states with mixed 



spin symmetry. For the moment, orbital relaxation due 
to KS self-consistency is ignored. The spin symmetries 
of these states can be explicitly determined through the 
action of the spin operators S 2 and S z , which yields (ig- 
noring orbital relaxations), 

S 2 |$ G ) = 0, S z |$ G ) = 



S 2 \<P T ) =2|$t>, 

S 2 |$Ml) = |*M1> 

S 2 |$M2) = 



^M2/ 



S z |$ T ) = 1 |*t) 

- |$M2>, S z |$Ml) = 

- |#mi), S z |$ M2 ) = (5) 



The above equations show that $mi,2 are riot eigen- 
vectors of the total spin operator, but rather a 
mixture of singlet and m z — triplet-states with 

(^mi^IS^I^mi^) = 1- Although not being eigcnfunc- 
tions of the actual many-body Hamiltonian, such mixed- 
states appear as variational extrcma in A-SCF calcula- 
tions since approximate exchange-correlation functionals 
depend on m z rather than on S 2 . Therefore, determin- 
ing the energ y o f the singlet excited state requires spin- 
purification |56l - [58| . whereby the energy of the singlet 
state is approximated through 



2E 



M 



(6) 



In this equation, Eg is the energy of the singlet state, 
Et is the energy of the triplet state calculated from the 
density corresponding to $t and -Em is the energy cor- 
responding to either of $mi or <!>M2 (note that they have 
the same energy, since the external potential is spin- 
independent). A brief derivation of Eq. ([5]) is provided 
in the Appendix. 

A commonly used approach for the calculation of the 
singlet excited state consists in performing the A-SCF 
calculation with the spin unpolarized (NSP) exchange- 
correlation functional [59l462{ . This type of calculation 
has shown to yield good agreement with experimental ex- 
citation energies. In the next section, we use and compare 
both the spin-purified approximation (SPA) and the NSP 
method for singlet energies. (In the Appendix, we pro- 
vide a brief discussion of the validity of both approaches.) 

As a consequence of the variational principle, the A- 
SCF calculation yields the lowest excited-state energy 
compatible with the imposed symmetry of the Slater de- 
terminant that was constructed from the one-electron 
KS states. Although it is recognized that the A-SCF 
approach lacks a formal justification (at variance with 
TDDFT that is based on the Runge-Gross theorem), the 
A-SCF technique is widely used in the literature and 
has been quite successful flTl - fl3| . One of the formal 
drawbacks of the A-SCF method is that the exchange- 
correlation functional is the same for both ground and 
excited states (i.e., the functional form is the same with 
the sole difference that the excited-state density is used 
in place of the ground-state density). More elaborate 
schemes with formal justification requires the excited- 
state calculation to be carried out by different function- 
als [HI, Hi| , by DFT methods with specific many-body 
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wavcfunction dependence [65|, [66| or by construction of 
ensembles with mixed ground and excited states [o7l - l69| . 
In either of these cases, the needed exchange-correlation 
functionals are unknown, and certain approximations are 
necessary. For constructing ensembles, one generally re- 
lies on statistical arguments extrapolated from the inter- 
acting electron gas [70]. The success of A-SCF for lowest 
energy excitations is plausibly related to the ground-state 
exchange-correlation functional being sufficiently accu- 
rate to describe excited-states in some systems. 



C. Hubbard model based functionals 

Typically, conventional approximate exchange- 
correlation functionals, such as PBE, yield a poor 
description of localized d states for transition metals. 
The DFT+U approach based on an orbital-dependent 
correction inspired from the Hubbard model has been 
widely used to obtain an accurate description of 
electronic localization in these cases p^H40l ]. More 
recently, the extension of DFT+U, to include inter-site 
interactions between atomic sites (DFT+U+V), has 
been quite successful, not only in systems with strong 
electronic localization (e.g. transition- metal oxides), but 
also in cases where electrons are delocalized (e.g., band 
insulators) [4l|. The V correction has also proved to 
be crucial in obtaining accurate structural properties 
of transition- metal dioxide molecules (7l| . and in de- 
scribing the dimerization of V atoms in VO2 accross the 
high-temperature to low-temperature phase transition 
through a DMFT approach (72|. In Ir complexes, while 
the on-site interaction U helps capturing the localization 
of d electrons on the central Ir atom, the inter-site V is 
expected to improve the description of the interactions 
with the surrounding organic ligands. 

In explicit form, the corrective functional in the 
DFT+U+V method is given by jU 

Euv = £^T3r[n"*(l-n'")] 

I, a 

-J2 ^Tr[„^„^] (7) 

I,J,a 

where U 1 and V 1 J arc on-site and inter-site interac- 
tion parameters respectively and indices / and J denote 
atomic sites. In Eq. (J7|), the occupation matrices n /Jcr 
are computed as 

= E fi tfmW) WlO (8) 
i 

where tp? are KS states, f[ are their occupations, and 
</>^j are atomic orbitals centered on site /. The first term 
in Eq. ([7]) that is proportional to U favors electronic 
localization on atomic sites, while the second term pro- 
portional to V leads to hybridization between orbitals 



on different sites, the ground-state configuration being 
eventually governed by the competition between the two 
opposite tendencies. 

IV. RESULTS AND DISCUSSIONS 

The structures of the three molecules [Ir(ppy)3, Fir- 
pie, and PQIr] have been optimized using the PBE func- 
tional, and kept fixed in all TDDFT and A-SCF calcu- 
lations. The TDDFT calculations are performed using 
PBE, B3LYP, and M06 functionals. 

The A-SCF calculations were performed using the 
GGA, GGA+U, and GGA+U+V functionals based upon 
PBE using the U and V parameters obtained from linear 
response [40]. The details of the linear- response calcu- 
lations for each of the molecules arc summarized in the 
following subsections. In some cases, re-optimizing the 
structure with forces due the U and V corrections was 
needed. 



A. Ir(ppy) 3 

The molecular structure of Ir(ppy)3 is shown in Fig. [2] 
The lowest triplet and singlet energies obtained within 




J 3 



FIG. 2: (Color online) Molecular structure of Ir(ppy)3 as a 
3d model and a chemical formula. 

TDDFT using the different semilocal and hybrid func- 
tionals are reported and compared with experiments in 
Table. HI Due to strong spin-orbit coupling present in 



TABLE I: Lowest TDDFT triplet and singlet excited state 
energies compared with experiment. All values are in eV and 
given relative to the ground-state energy. 





Exp 


GGA 


B3LYP 


M06 


T 


2.4" 


2.09 


2.55 


2.55 


S 


2.6 6 -2.7 c 


2.16 


2.75 


2.79 



"From Refs. 0, I73l482l , from the highest peak in the phosphores- 
cence spectra. 
6 From Ref. [1 

c From Refs. \tA |78|. |82| from the second peak/shoulder in the ab- 
sorption spectrum (the first absorption feature is assumed to be 
due to the triplet state). 
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these systems, it is in principle not possible to distin- 
guish between singlet and triplet states since the actual 
cigenstatcs are mixtures of them. However, the triplet 
energy can be measured quite accurately in experiment 
from the phosphorescence spectra. The largest peak in 
the phosphorescence spectrum determines the energy of 
the triplet state, and can be determined with high accu- 
racy [3l. 1 73148 2j . Nevertheless, in contrast to the triplet 
state, the experimental determination of the singlet en- 
ergy is very difficult. One could in principle measure 
the onset of optical absorption of singlet excitons for this 
purpose. However, due to strong spin-orbit coupling, the 
optical absorption starts at much lower energies than the 
actual singlet energy, leaving a long tail that extends 
deeply into the high- wavelength region (77T483T ]. An- 
other possibility is to determine the strongest peaks in 
the absorption spectra, which would yield much higher 
singlet-state energies than the lowest one. Alternatively, 
the LUMO-HOMO gap when the LUMO is measured by 
optical absorption (the so-called optical LUMO) has also 
been used to determine the singlet state energy [§3, [85[ . 
It is important to note that the optical LUMO contains 
the binding energy of the exciton unlike an inverse photo- 
emission measurement where the exciton binding energy 
is eliminated [8(1 [8?]] ■ However, the strong spin-orbit cou- 
pling also affects such measurements since they also rely 
on optical absorption. Therefore, comparison of the cal- 
culated lowest singlet energy with the experimental data 
is not always valid, and one needs in principle to compute 
the absorption spectrum including spin-orbit coupling ef- 
fects, which can be computationally very expensive. Rel- 
ativistic effects have only recently been studied for similar 
dyes and was shown to result in better agreement with 
experiment [88j . Another uncertainty that might affect 
the analysis of the experimental data is whether the mea- 
surements are performed in the gas phase, in solution or 
for the molecule absorbed on the surface of a thin film. 
However, the correspondence between different types of 
measurements are well-known in the literature [891 ] . 

In spite of recognized difficulties in making comparison 
with optical measurements, we report the calculated low- 
est singlet excited-state energies (and trends in splittings 
from the triplet excited states) because they are critical 
for the purpose of guiding the design of OSCs with Ir 
dyes used as sensitizers Q. 

As can be seen from Table Q] the triplet excited- 
state energy is predicted to be within 0.1 eV from the 
experimental value by the hybrid functionals (B3LYP 
and M06), whereas it is notably underestimated by the 
semilocal functional (PBE). Although as already men- 
tioned, it is very difficult to determine the lowest singlet 
energy experimentally, the hybrid functionals yield pre- 
dictions within an error of 0.2 eV relative to the reported 
experimental values, whereas PBE again underestimates 
it. The excited-state energies reported in Table fl] are also 
in^agre ement with the previous results in the literature 

In order to gain further insight into computational pre- 



dictions, we depict the probability density of the HOMO 
and the LUMO for the ground state within PBE in Fig. [3] 
As can be seen, the HOMO is predominantly a metal- 




FIG. 3: (Color online) Ir(ppy) 3 HOMO(left) and 
LUMO (right) calculated using the ground-state electronic 
density. 

centered state with some contribution from the ligands 
while the LUMO is almost entirely ligand-centered, in- 
dicating that optical excitations correlate with metal-to- 
ligand charge-transfer (MLCT) processes, as was previ- 
ously recognized in the literature [90l - l95| . This charge- 
transfer property of the excited states also explains why 
the local and semilocal adiabatic approximation under- 
estimate their TDDFT energies. 

The calculation of U and V parameters requires per- 
turbing each atomic site separately for the construction 
of the response matrices. Instead of performing such a 
computationally demanding calculation including the re- 
sponse of all incquivalent sites, we concentrate only on 
the central Ir atom and its nearest neighbor C and N 
atoms. This approach is justified by our expectation that 
the electronic correlations are only important for elec- 
trons localized on the Ir atom. The atoms further away 
from Ir in the organic ligands (i.e., the distant C and H 
atoms) are taken into account in an average manner as 
a charge reservoir centered at the Ir site, as explained in 
Ref.[96|. The linear response calculation yields the on-site 
U parameter for Iv-d (Ud) and Ir-s (U s ) states, the V pa- 
rameter between Ir-d,s and neighboring C-p,s and N-p,s 
states [V(lrd, s , N p _ s ) and V(lrd, s ,C PyS )] and the on-site 
interaction parameter between Ir d and s states (Uon)- In 
Table. Ullwe report only a subset of these values, which af- 
fects the electronic structure most significantly. Further 
calculations, not reported here, have proven that other 
interaction parameters have no effect on the excited- 
state energies. The reported values of V^Ir^, N p ^ s ) and 
V(lrd, Cp,s) ar c average values, since the distances be- 
tween the three N atoms and the three C atoms to Ir 
atom differ slightly from one ligand to another. How- 
ever, the differences are within the numerical precision 
of the linear-response calculations. We have also com- 
puted the inter-site interaction parameters between the 
Ir-d states and C-p, s states which are not nearest neigh- 
bors to the Ir atom. This calculation is performed by 
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TABLE II: Linear-response values of U and V parameters. 
We use the notation V(Ir d , N p>s ) = VjL PtS and V^Ir^, C p , s ) = 
V d c _ v s . All values are in eV. 



u d 


V N 


V? 

v d—p 


V d N - s 


v£L a 


7.36 


1.62 


1.32 


2.39 


2.65 



isolating a chain of atoms which starts from the Ir atom 
and ends at a C atom in one of the ligands. The rest of 
the atoms in the molecule is treated as a charge reservoir, 
which is placed at the Ir site. As can be seen from the 
results shown in Figs. 01 and [SI the inter-site V parame- 
ters decay as a function of the distance from the Ir atom. 
This result provides a quantitative justification for the 
use of the nearest neighbor V parameters between Ir-N 
and Ir-C pairs only. 



> 

> 




FIG. 4: Inter-site V parameters calculated for a Ir-C-C-C-C 
chain, as a function of distance. V dp denotes the interaction 
parameter between Ir d and C p states, while V da denotes that 
between Ir d and C s states. 



> 




FIG. 5: Inter-site V parameters calculated for a Ir-N-C-C- 
C chain, as a function of distance. V dp denote interaction 
parameter between Ir d and C,N p states, while V ds denote 
those between Ir d and C.N s states. 



The results of the A-SCF calculations using both the 
SPA singlet-state energy in Eq. and its NSP coun- 
terpart are reported in Tabic Mil Compared with the 



TABLE III: A-SCF calculation of the lowest triplet and sin- 
glet states for Ir(ppy)3. All values are in eV and measured 
from the ground-state energy. 





GGA 


U 


u+v 


(U+V) rcl 


Exp 


T 


2.32 


2.27 


2.41 


2.44 


2.4 


M 


2.34 


2.29 


2.43 


2.46 




S [NSP] 


2.50 


2.90 


2.78 


2.73 


2.6-2.7 


AE TS [SPA] 


0.05 


0.04 


0.02 


0.04 




AE TS [NSP] 


0.18 


0.62 


0.37 


0.29 


0.2-0.3 



TDDFT calculations reported in Table U the triplet en- 
ergy resulting from the GGA+U and GGA+U+V cal- 
culations are slightly smaller than the triplet state en- 
ergy calculated with hybrid functionals. Notice that the 
GGA functional within the A-SCF approach yields re- 
sults closer to experiments than TDDFT(PBE). This is 
due to the ability of A-SCF to partially rectify charge- 
transfer errors. Although the energy of the triplet-state 
was obtained accurately, when the singlet energy is eval- 
uated through the spin-purification formula of Eq. ([5]), 
the A-SCF singlct-triplct splitting is almost vanishing. 
On the other hand, the NSP calculation of the sin- 
glet state (S) yields results in better agreement with 
experiment and TDDFT. In quantitative terms within 
NSP, GGA+U overestimates the singlet energy by 0.3 eV, 
whereas GGA+U+V improves the agreement to an ac- 
curacy of approximately 0.2 eV. This improvement re- 
sults from the inclusion of inter-site interactions through 
V, which corrects the overlocalization resulting from the 
straight use of the on-site Ud- In addition, when the 
molecular structure is re-optimized with the GGA+U+V 
functional [i.e. with additional forces coming from the 
corrective terms in Eq. (J7J] , the triplet-state almost 
matches with the experimental value, while the singlet 
energy is overestimated by only 0.1 eV. Although both 
TDDFT and the Hubbard-corrected functionals yield re- 
sults that are close to experiments, the experimental val- 
ues themselves have an accuracy close to 0.2 eV (85[. Er- 
rors in the singlet energies are in principle much larger 
due to problems inherent to their experimental determi- 
nation, as discussed above. 



B. Fir pic 

The molecular structure of FIrpic is shown in Fig. [6] 
The lowest triplet and singlet excited-state energies cal- 
culated with TDDFT using PBE, B3LYP, and M06 func- 
tionals are also compared to experimental data in Ta- 
ble IIVI A salient feature in this comparison is the fact 
that the PBE functional underestimates the energies of 
both excited states similarly to Ir(ppy)3 (cf. Table UJ. 
Instead, hybrid B3LYP and M06 show a significant im- 
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FIG. 6: (Color online) Molecular structure of FIrpic as a 3d 
model and a chemical formula. 



and V interaction parameters were calculated using the 
approach described in subsection IIV Al We have also 
checked the validity of this approach by isolating several 
chains containing Ir, C, N and O atoms and computing 
V between pairs that are not nearest neighbors. We ver- 
ified that the interaction parameters between the Ir and 
its nearest neighbor atoms are important, while the in- 
teraction parameters between Ir and more distant shells 
of neighbors vanish similarly as in Figs. [4] and [5] Due to 
this strong similarity of the results with Ir(ppy)3, we do 
not report these calculations here. 

The calculated U and V interaction parameters for 
FIrpic are reported in Table. |Vj We report only the 



TABLE IV: Lowest triplet and singlet excited-state energies 
obtained with TDDFT and from experiments for FIrpic. All 
values are in eV and measured from the ground-state energy. 



TABLE V: Linear-response values of U and V parameters. 
All values are in eV. 





Exp 


GGA 


B3LYP 


M06 


T 


2.6° 


2.14 


2.67 


2.66 


S 


3.3 6 -2.9 c 


2.24 


2.94 


3.00 



"From Refs. 0, 197l4l02l from the highest peak in the phosphores- 
cence spectrum. 

'From Ref. l97l, from the second peak in the absorption spectrum. 
'From Ref. Il02l , from the onset of optical absorption. 



provement, and yield a triplet energy very close to ex- 
perimental value. Despite the difficulties discussed in the 
previous section, the hybrid functionals are always within 
the experimental uncertainty. Therefore, we argue that 
they are fairly accurate. 

To further explore the properties of the excited states, 
the HOMO and the LUMO densities in the ground state 
of FIrpic, are plotted in Fig. [7] within PBE. Similarly to 





FIG. 7: (Color online) FIrpic HOMO(left) and LUMO(right) 
calculated using the ground-state electronic density. 

Ir(ppy)3, the HOMO is mostly located on Ir with consid- 
erable contribution from the ligands, whereas the LUMO 
is almost entirely ligand-centered on the attached chro- 
mophores [predominantly the picolinate (pic) ligand], 
with some contribution coming from the metal center. 

The A-SCF calculations were performed using the 
GGA, GGA+U, and GGA+U+V functionals. The U 



u d 


v d—p 


v d — p 


V? 

' d — p 


V d N -s 






7.17 


1.69 


1.29 


2.0 


3.06 


2.56 


5.77 



parameters that affect the electronic structure of the 
molecule most significantly, and disregard the ones that 
do not contribute at all. Notice that the parameters Ud, 
U(Ir d , C p ), and V(lv d , N p ) are very close to the ones 
which are calculated for Ir(ppy)3 reported in Table. HH 

The results of the A-SCF calculations using both the 
SPA and NSP calculations for the singlet state are re- 
ported in Table EH As for Ir(ppy) 3 , GGA within A-SCF 



TABLE VI: A-SCF calculation of the lowest triplet and sin- 
glet states. All values are in eV and given relative to the 
ground-state. 



GGA 



U 



U+V (U+V) rc i Exp 



T 

M 
S [NSP] 
AE TS [SPA] 
AE TS [NSP] 



2.46 
2.50 
2.65 
0.08 
0.18 



2.79 
N.A. a 

3.00 
N.A" 

0.20 



2.56 
2.60 
2.88 
0.08 
0.31 



2.52 
2.56 
2.81 
0.08 
0.29 



2.6 
3.3-2.9 
0.3-0.7 



"No convergence is achieved for these calculations. 

yields excited-state energies in better agreement with 
experimental data compared to GGA within TDDFT. 
While GGA+U overestimates the triplet energy by ap- 
proximately 0.2 eV, GGA+U+V lowers it by 0.2 eV, that 
is, within 0.1 eV of experiments. Additionally, the singlet 
state energies predicted by GGA+U and GGA+U+V 
fall within the experimental range, as was the case for 
TDDFT with hybrid functionals. 

Yet, compared to Ir(ppy)3, the excited-state energies 
are higher in FIrpic, explaining the blue phosphorescence 
of FIrpic. We thus confirm a well-known experimental 
result that the attachment of F atoms in the ligands sta- 
bilize the HOMO level and increase the LUMO-HOMO 
energy gap [l03lll04l |. Since F is highly electronegative, 
it decreases the amount of charge localized on the metal 
center, thereby reducing the total Coulomb repulsion felt 
by these electrons and stabilizing the HOMO level. This 
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comparison with Ir(ppy)3 suggests that further adjust- 
ments on the triplet emission energy can be obtained 
by functionalizing the (pic) ligand, since the LUMO is 
mainly localized on it, as can be seen from Fig. [7J 

C. PQIr 

The molecular structure of PQIr is shown in Fig. M 
The lowest triplet and singlet excited-state energies 




FIG. 8: (Color online) Molecular structure of PQIr as a 3d 
model and a chemical formula. 

within TDDFT using GGA (PBE), B3LYP, and M06 are 
reported in Table. IVIII As for the other molecules stud- 



TABLE VII: Lowest triplet and singlet excited-state energies 
obtained with TDDFT and from experiments for PQIr. All 
values are in eV and measured from the ground state. 





Exp 


GGA 


B3LYP 


M06 


T 


2.1 


1.63 


2.02 


2.01 


S 


2.3 6 


1.74 


2.35 


2.37 



"From Refs. [^, |75| . Il05l from the highest peak in the phosphores- 
cence spectr um. 

''From Ref. [l05l4l07l. from the LUMO-HOMO energy difference, 
where the LUMO is treated as an "optical LUMO" which also 
contains the exciton binding energy [84l . |85|| (which is a singlet 
state) . 

ied in this paper, TDDFT(PBE) underestimates both the 
triplet and the singlet energies while the hybrid func- 
tionals yield results within an error of 0.1 eV relative to 
experiment. 

The HOMO and the LUMO charge distributions are 
depicted in Fig.[9]for the ground state within PBE. Here, 
in accordance with the other Ir dyes, the HOMO is pre- 
dominantly centered on the metal, with some contribu- 
tions from the ligands. Instead, the LUMO is almost en- 
tirely ligand-centered. Notice that the tetramcthyl hep- 
tancdionate (tmd) ligand only contributes (marginally) 
to the HOMO, and the LUMO is concentrated on the 
larger phenyl-quinoline (pq) ligand. 

The calculated U and V parameters are reported in 
Table IVIIIl As before, we only report interaction pa- 
rameters that affect the electronic structure most sig- 
nificantly. The validity of restricting our calculations to 




FIG. 9: (Color online) PQIr HOMO(left) and LUMO(right) 
calculated using the ground-state electronic density. 



TABLE VIII: Linear-response values of U and V parameters. 
All values are in eV. 



u d 


v d—p 


v c 


v d — p 


Vts 






7.26 


1.71 


1.39 


1.79 


2.66 


4.43 


5.40 



Ir and its nearest neighbors was also verified for PQIr, 
by computing U and V parameters between the atoms 
along several "chains" from the central Ir towards the 
external ligands and showing that only nearest neighbor 
interactions are important. The calculated Ud and V for 
Ir-N, Ir-C, and Ir-0 pairs between d and p states are also 
very close to the ones computed previously for the other 
two molecules. This finding suggests that these parame- 
ters could be considered to be universal and used for fu- 
ture studies without recomputing them for each molecule 
(provided that the same pseudo-potential and exchange- 
correlation functional are used). 

The results of the A-SCF calculations using the GGA, 
GGA+U, and GGA+U+V functionals are reported in 
Table Hxl Again, GGA yields results in better accordance 



TABLE IX: A-SCF calculation of the lowest triplet and sin- 
glet states for PQIr. All values are in eV and given relative 
to the ground-state. 





GGA 


U 


u+v 


(U+V) rol 


Exp 


T 


1.81 


2.11 


1.93 


1.90 


2.1 


M 


1.79 


2.18 


1.98 


1.95 




S [NSP] 


2.01 


2.34 


2.23 


2.18 


2.3 


AE TS [SPA] 


0.06 


0.14 


0.10 


0.10 




AEts [NSP] 


0.19 


0.23 


0.30 


0.28 


0.2 



with experiment within A-SCF than within TDDFT. 
The best agreement with the reported experimental val- 
ues in Table IVIII are obtained with GGA+U, using the 
NSP calculation of the singlet state. This might seem to 
contradict the results for Ir(ppy) 3 and FIrpic, where the 
best agreement was found with the inclusion of V. How- 
ever, it should be noted that the experimental results 
have an accuracy of approximately 0.2 eV [H[. Thus, 
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GGA+U+V results are within this range of this accu- 
racy. 

The excited-state energies of PQIr are smaller than 
those of Ir(ppy)3 and FIrpic, which accounts for its red 
phosphorescent emission. This is due to the stabiliza- 
tion of the LUMO level resulting from the larger size 
of the (pq) ligand where it is mostly centered on (see 
Fig. In fact the (pq) ligand has an additional ben- 
zene ring compared to Ir(ppy) 3 . A larger ligand structure 
can be expected to lower the electron-electron interac- 
tions and reduce the electronic kinetic energy through a 
more pronounced derealization, hence the stabilization 
of the elec tronic states it hosts (the LUMO level in this 
case) [tol l. 



D. Comparison of TDDFT and A-SCF 

As discussed briefly in subscction llll A\ TDDFT yields 
inaccurate energies for charge-transfer excitations, specif- 
ically if used with local and semilocal functionals, like 
GGA(PBE). As we have shown in Figs. C3 and H the 
excited states of the Ir dyes studied in this work are 
of MLCT type. The MLCT character of the excited- 
states explains why GGA underestimates the excited- 
state energies, whereas the hybrid functionals provide 
better agreement with the experimental measurements, 
which are reported in Tables HI HV| and IVHI Instead, 
PBE within the A-SCF approach using the NSP repre- 
sentation of the singlet state proves to be more accu- 
rate than TDDFT, as can be seen from Tables Hill ED 
and IIX1 Indeed, the fact that A-SCF could perform 
better than TDDFT for charge-transfer and Rydberg 
excitations is already known from the literature [Tip 
[H, [H HU- TDDFT relies on KS energies determined 
from the ground-state configuration, and in the case of 
charge-transfer excitations, the effect of unoccupied or- 
bitals are not accurately taken into account by semilo- 
cal functionals since the off-diagonal block matrices (M) 
nearly vanishes in Eq. (TIJ. This is due to the fact that 
the exchange-correlation potential decays exponentially 
with distance r, rather than as 0(r _1 ). Instead, A-SCF 
takes into account the orbital relaxation of ground-state 
KS states, since those are reconstructed at each iteration 
from an excited-state density [Eq. ((4])] . If expanded us- 
ing the basis of unrelaxed orbitals, the relaxed orbitals 
include components from the unoccupied manifolds that 
are not properly taken into account by TDDFT used with 
(semi)local functionals. As a result, the effect of the un- 
occupied manifold (with respect to ground-state KS or- 
bitals) is partially included in A-SCF. Along the same 
lines, Ref. (5J has shown that the long-range behavior of 
the A-SCF energies can be more accurate than TDDFT. 

Comparing Tables U El EED with Tables UTTl [VTl HXl 
one can observe that hybrid functionals within TDDFT 
performs similarly to A-SCF calculations with Hubbard 
corrections. This shows that both methods can be used 
interchangeably In explicit terms, one can use A-SCF 



with hybrid functionals or TDDFT with Hubbard U and 

V corrections, where the latter approach would allow for 
significantly less computational effort for large systems 
than hybrid functionals. One other notable advantage 
of the Hubbard corrected functionals is that the U and 

V parameters directly control the amount of electronic 
localization on the Ir atom and the hybridization be- 
tween Ir and the surrounding ligands. This knowledge 
could provide a functionalization strategy for tuning the 
excited-state energy levels, as discussed briefly in the 
next subsection. 



E. Effects of U and V 

In order to investigate the effects of the Hubbard cor- 
rections more systematically, we provide a simple model 
of the system, based on two atomic states, one centered 
on the metal and the other centered on the ligand, as 
shown schematically in Fig. [TDJ In this simplified pic- 




FIG. 10: A simple description of the MLCT character of the 
excited states in the Ir complexes. 

ture, we assume that the HOMO levels i/jJj have large 
overlap with the metal orbital (f> M , and a small overlap 
with the ligand orbital <fi L . Instead, the LUMO levels 
have large overlap with the ligand orbital <p L but a small 
overlap with the metal orbital <f> . The contribution of 
the Hubbard corrective energy in Eq. (J7J) in the case 
of two atomic orbitals cf> M and (j> L (which are assumed 
non-degenerate) becomes 

w Um \~* Ma Ma\ irML MLo LMa 

i/uv = — — 2_^i n (J- — n J — V n n 

a a 

(9) 

where 

n M ° = Ylfi\(4> M m)\ 2 

i 

n ML* = Y,i:^ M mm4> L ) (io) 

i 

In the above equations, the occupation numbers f? are 
either 1 or 0. We also assume that the occupation ma- 
trices are symmetric, i.e., n MLa = n LMa ' . We can con- 
struct the occupation matrices n M , n L , n ML by ignoring 
the orbital relaxation and using the levels schematically 
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represented in Fig. [TJ For example, in the triplet state, 
the occupation matrices are given by 

Aft Aft , A t ATI Ml A I 

n T = n G + A n ' N+1 , n T * = n G * — A n N 

MLt _ ML\ a ML\ MLl _ MLl _ * MLl 

i irp — G i 7V-|-1 ' T — G N 

(ii) 

where 

A< +1 = |^ A/ |Vjv + i)l 2 > Ani^K^I^)! 2 

An££ ss (^|^ +1 ><VUil^> 

An^ = (0 M |^)(^I^> (12) 

Since "0^ is dominantly metal-centered, while V'jv+i is lig- 
and centered, An] v+1 1 is satisfied, whereas An N+ [ 
and An^ /L ^ can in principle be larger, but still consider- 
ably smaller than 1. Instead, An^ is of the order of 1. In 
Eq. lfTTj) . the subscripts G and T refers to occupation ma- 
trices evaluated from the ground-state and triplet-state, 
respectively. We can calculate the contribution to the 
excited-state energies from Hubbard terms as E]j V — E^ v 
by using Eqs. ([§])- (fT2"]) . The contribution coming from 
the on-site interactions [i.e., the first term in Eq. ©] is 
given by 

E%-E®=ui> = Ani (l + Ani - rg) 

+ U -f An\ +X (1 - n£) + O (An% +1 ) 

(13) 

where we have ignored quadratic terms in Anj v+1 and 
used the fact that n G = n G = \n G . The first term in 

Eq. (JT3J) is negative, since 1 + A l N > 1, while n£f < 1. 
Instead, the second term is positive. Therefore, Eq. (fT"3|) 
shows that the excited-state energy decreases with in- 
creasing Um, which is expected since larger Um destabi- 
lizes the HOMO level by adding larger penalty when it is 
doubly occupied in the ground-state. However, for cases 
where there is considerable metal and ligand overlap, the 
second term in Eq. (fT3"|) can reverse this effect through 
a larger Anj v+1 . To test this prediction, we have plot- 
ted in Fig. [TT] the energies of the triplet and singlet state 
(obtained from a NSP calculation) as a function of Ud for 
Ir(ppy)3. As can be seen from Fig. [ITJ the triplet energy 
increases for Um < 5 cV, but decreases for larger values 
of Um- This decreasing behavior is correctly predicted 
by Eq. (|13[) . Instead, the increase in the triplet energy 
for Um < 5 eV is probably a result of important modifi- 
cations in the KS states due to orbital relaxation which 
is implicitly ignored in deriving Eq. (|13|) . For the case 
of the molecules studied in this work, Ud ^ 7 eV. There- 
fore an increasing Ud would decrease the triplet energy, 
mainly due to the destabilization of the HOMO level. 
The singlet state is a result of a NSP calculation, and its 
behavior as a function of Um is thus difficult to assess 




23456789 
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FIG. 11: Triplet and singlet energies computed as a function 
of on-site U on Ir d states. No inter-site V is included. 

since in principle it corresponds to an ensemble which 
contains an arbitrary combination of possible different 
Slater determinants, as discussed in the Appendix. 

The contribution from the inter-site interaction terms 
in Eq. (|9]) to the excited-state energy can be calculated 
similarly. Doing so, we obtain 

E%-E$ = u% = -V ML n% L (An^I-An^) 

+0 (A 2 ) (14) 

where we have ignored the quadratic terms in Anjv and 
Arijv+i. Note that the term given in Eq. (|14p could 
be either positive or negative depending on the relative 
magnitude of the projections An^f_[ and An^ L ^ . For 
example, in the case of a strong metal-LUMO overlap but 
a smaller ligand-HOMO overlap, V decreases the triplet 
energy. Instead, in the opposite case where the ligand- 
HOMO overlap is larger than the metal-LUMO overlap, 
V increases the triplet energy. In order to validate these 
predictions, we have plotted the triplet and singlet (from 
NSP calculation) state energies as a function of the V 
between Ir d and C and N p and s states in Fig. Q21 In 
these calculations, we have fixed the value of Ud to the 
calculated value of 7.36 eV and used the same V between 
all the Ir-C and Ir-N pair interactions. As can be seen 
from Fig. |T2l the triplet energy increases for V < 1.5 eV 
and decreases for V > 1.5 eV, reflecting both type of be- 
haviors shown in Eq. (fT4|) . In this case, the orbital relax- 
ation effects, which are ignored in obtaining Eq. (|14p are 
more critical since they affect the overlap of the HOMO 
and LUMO levels (V>jy, V'at+i) w i tn tne metal-and-ligand- 
centered states. This contribution can very easily invert 
the behavior of the triplet energy as a function of V. In 
summary, Figs. [TT] and [12] show that the excited-state 
energies critically depend on the U and V values, which 
highlights the sensitivity of these energies to small chem- 
ical changes. In fact, as discussed in previous sections, a 
change in the chemical composition of the ligand groups 
can result in modifications to the electronic structure of 
these molecules analogous to those obtained by varying 
U and V. This conclusion was also obtained in a different 
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V(lr d ,NC p , s ) [eV] 

FIG. 12: Triplet and singlet energies computed as a function 
of inter-site V between Ir d and C,N s,p states. U = 7.36 eV 
fixed at the computed value. 

way by using model Hamiltonians in Refs. 1 1091 and IllOi 
In fact, the results presented in Tables [TTl IVl and fVIIII can 
be used as effective parameters in a modcl-Hamiltonian 
study for these systems. 

Summarizing, for the three molecules that we have 
studied in this work, two distinct behaviors are observed. 
For Ir(ppy) 3 , the effect of the U and V corrections is 
mainly to tune the HOMO level. The addition of U de- 
creases the excited-state energies (see Table IIII[) due to 
the destabilization of the HOMO level, which is domi- 
nantly metal-centered. This destabilization arises from 
an increased Coulomb repulsion between the two elec- 
trons occupying the HOMO level. Instead, with the addi- 
tion of V, the HOMO level is stabilized, increasing the ex- 
citation energies. This is due to the increased hybridiza- 
tion between Ir and neighboring atoms, which depletes 
the Ir d states, lowering the Coulomb repulsion on them. 
In contrast, for FIrpic and PQIr, the addition of U and 
V affects mainly the LUMO level. Including U increases 
the excited-state energies (see Tables |VI] and |FX|) due to 
the destabilization of the LUMO level. This trend can 
be ascribed to the fact that the LUMO level contains a 
non-negligible contribution from the metal, as compared 
to Ir(ppy) 3 (compare Fig. [3] with Figs. [7] and [9]). Instead, 
the addition of V stabilizes the LUMO by lowering the 
fraction of electrons localized on the metal-center. In 
a nutshell, this study suggests that the substitution of 
atoms close to the Ir center can be as effective in the 
functionalization of the complex as the well known addi- 
tion of F or other electronegative species to the external 
part of the ligands or modifications of their structure. 



namely, GGA+U and GGA+U+V. The results obtained 
in both approaches are in good agreement with experi- 
ment. The Hubbard corrections U and V were computed 
from ab-initio calculations and provide a measure for lo- 
calization and hybridization of Ir d states with neighbor- 
ing organic ligands in the ground and excited-states. This 
knowledge is used to infer possible strategies for tuning 
the excited-state energies of the studied molecules. The 
gained insight underscores the interest of Hubbard cor- 
rections in unveiling the electronic origin of dye phos- 
phorescence. In addition, we have also investigated the 
validity of the spin purification (SPA) and the nonspin- 
polarized (NSP) calculations for the computation of the 
singlet excited-state energy. We found that the SPA ap- 
proach (based on the construction of a state with mixed 
spin symmetry) clearly underestimates the experimen- 
tal singlet energies. Instead, the NSP calculation yields 
much more accurate results. While the failure of the spin 
purification formula can be understood from the inad- 
equacies of the conventional exchange-correlation func- 
tionals, the remarkable success of the NSP approach in 
capturing the singlet remains a relevant open question. 

The present work can be considered as the starting 
point of two main research directions involving method- 
ological developments on one side and the design and op- 
timization of better molecular complexes, on the other. 
From a methodological point of view this study high- 
lights the necessity to improve functionals that can dis- 
tinguish between different (S 2 ) configurations (with the 
inclusion of spin-orbit interactions). On the application 
side this work establishes DFT+U+V as a valid compu- 
tational tool (thanks to its ability to capture the charge- 
transfer character of excited states) to efficiently screen 
useful modification to the ligand structure of the con- 
sidered molecules that could improve their performance. 
Furthermore, the dependence of the energy splittings on 
the values of U and, more importantly, of V, suggests 
that the choice of the chemical species directly bonded 
with the Ir center and of the structure of the ligands in its 
coordination shell may represent another valuable route 
to functionalization that can be explored in addition to 
the modification of their most external part. 
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Appendix A: Comparison of spin-purification and 
NSP calculations 



In section ITVl we have shown that the A-SCF energy 
for the mixed spin states 3>mi,2 used in the SPA formula 
Eq. yields much lower singlet state energies than the 
experimental results. Instead, the NSP A-SCF calcula- 
tion yields results within the experimental accuracy of 
0.2 eV. In this Appendix, we provide some insight into 
this finding. 

Approximate exchange-correlation functionals depend 
only on the spin component m z , and cannot distinguish 
between m z = triplet and singlet states. Moreover, 
with approximate exchange-correlation functionals and 
integer occupations, only configurations that can be rep- 
resented by a single Slater determinant are obtained as 
stationary densities in A-SCF calculations. Instead, the 
m z = triplet and singlet states correspond to multi- 
Slater-detcrminant configurations, represented by the fol- 
lowing density matrices (excluding orbital relaxation ef- 
fects) 



|$Ml)($Ml| - |*Ml)($M2| 
-|$M2>($Ml| + |$M2)($M2| 



(Al) 



A>,0 



1 r 

2 



|$Ml)($Ml + |$Ml)($M2| 
+ |$M2>($Ml| + |$M2)($M2| (A2) 



where Di corresponds to the density matrix of the m z = 
triplet, while -Do,o corresponds to the density matrix 
of the singlet state. Eqs. (|A1[) and (|A2[) can easily be 
verified by evaluating the expectation value of the total 
spin square operator, and making use of Eq. (J5]) as 



(S ?2 ) 1 ,o = Tr 




a V = T* D 0fi S 2 =0. 



(A3) 



Instead, the single-particle density is the same for both 
states, and is given by 

n(r) = Tr Dh(r 

= ^ ($Ml|"(r)|$Ml) + r ( $ M2Kr)|$M2)- 

(A4) 

Notice that the off-diagonal terms ($mi|"-(i')|<I > M2) van- 
ish when orbital relaxation effects arc ignored (i.e. the 
states $mi,2 are fixed). With the inclusion of orbital 
relaxation effects, the densities of m z = triplet and 
singlet states would be different. However, realization 
of this difference requires an ensemble dependent func- 
tional (67j-f70| . More specifically, the ensemble exchange- 
correlation functional should depend on the weights of 
the Slater determinants appearing in the density matri- 
ces Eq. (|Aip and (|A2[) . In fact, any density matrix of 



the form qbi$ + (1 — q) Dq.o with q < 1, results pre- 
cisely in the same one particle density Eq. (|A4|) . with 
(S 2 ) — 2q. Without a functional that depends on the 
ensemble weights (i.e., q in this case) or a functional that 
can distinguish between states with different S 2 , it is 
not possible to know which ensemble n(r) corresponds 
to. In the absence of such a functional, one relies on 
the calculation of the expectation value of the many- 
body Hamiltonian operator using DFT with standard 
exchange-correlation functionals, which are given by 



E lfi =Tr D li0 U 
1 



2 ($Ml|ft|$Ml) + ^ (*M2|H|$M2) 



1 



Ea.o = Tr 



[($Ml|H|$M 2 ) +C.C.] 



Dq o % 



- ($Ml|H|$Ml) + ^ ( < J , M2|H|$M2) 
+ \ [(*Ml|W|*M2) + C.C.] 



(A5) 



Assuming that the expectation value of the many-body 
Hamiltonian can approximately be given by A-SCF cal- 
culations, one can identify the mixed-state energy as the 
variational extrema given by 



(*M1,2|^|*M1,2> = E M . 



(A6) 



Since the Hamiltonian does not contain any external 
magnetic fields, m z = 1 and m z = triplet states should 
have the same energy, which requires 



($t|"H|$t) = E T = Ei. a 



Then, adding the two equations in Eq 
off-diagonal terms, and using Eqs 
obtains the spin purification formula 
E T + E S = 2 Em (with E s = E ,o) K 



(A7) 



5|) to cancel the 
7J) and ([A6]). one 
iven in Eq. ([5]): 
5? 



The reason why the A-SCF energy of the mixed-state 
(.Em), used in spin-purification formula Eq. ([5]), under- 
estimates the experimental singlet energies can be un- 
derstood, in part, by carefully examining the assump- 
tions that were made in the discussion above. First, the 
assumption in Eq. (|A7[) is not justified. Approximate 
ground-state exchange-correlation functionals can distin- 
guish between states with different m z , but not with dif- 
ferent (S 2 ), as discussed previously. Therefore, Eq. (|A7[) 
assumes the equality of energies of a single determinant 
state $t and an ensemble £>i,o, which have different m z . 
In other words, ($t|%|$t) 1S n °t a g°°d estimate of 
Et. Such an identification is clearly beyond the capa- 
bilities of approximate exchange-correlation functionals. 
Eq. (| A7|) can be justified only if an appropriate ensem- 
ble exchange-correlation functional is used. Second, the 
spin purification formula is plagued by the negligence of 
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orbital relaxation effects. The ensemble density matrices 
of the singlet and m z = triplet states of Eqs. (|A1[) and 
(|A2j) suggests that they are equal- weight combinations of 
the mixed-states $mi and <I>M2- The spin purification for- 
mula critically depends on this equivalence of the weights. 
However, when orbital relaxation effects are included, KS 
orbitals corresponding to $mi arc $M2 are different and, 
in general, the transformation between them requires an 
infinite expansion, including infinite set of states in the 
unoccupied manifolds. Thus, an actual ensemble repre- 
sentation of the singlet and m z = triplet-states should 
contain contribution from such an infinite expansion (e.g. 
on unrelaxed KS orbitals). Namely, the ensembles must 
include a linear combination of many excited-state de- 
terminants with orbitals ip% occupied where a > N + 1. 
Such an approach is not viable, since constructing a func- 
tional that depends on the weights of each excited Slater 
determinant is not easy. At this point we would like to 
stress that ensemble DFT of Refs. l67l - [70l does not suf- 
fer from such problems, provided that an ensemble de- 
pendent exchange-correlation functional exists. In this 
case, the energy of each state in the ensemble is uniquely 
determined by the Hohenberg-Kohn theorem for ensem- 
bles (67|. 

The success of the NSP A-SCF calculation is more 
difficult to understand. The NSP state used for the A- 
SCF calculation, $nsPj exactly corresponds to a transi- 
tion state with orbitals ipN and ipN+i occupied by half 
spin-up and half-spin down electrons, as shown schemat- 
ically in Fig. 1131 It is well-known that the single parti- 




NSP 



1/2 1/2 

t ¥ 

f Nh 

1/2 1/2 

■v., 



3>NSP 



FIG. 13: "^nsp schematically represented both by spin- 
unpolarized orbitals and an equivalent transition state. 

cle density of a transit ion s tate [lll| is equivalent to an 
ensemble density [H, Ill2| . Indeed, the corresponding 
single particle density obtained from "1>nsp is given by 



iV-l 



n NSP (r) = ICO) 



l4r(r)l 



+ |^(r)| 2 + |^ + iW| 2 + l^ + i(r)| 2 

(A8) 

When orbital relaxation effects are taken into account, 
Eq. (|A8|) still holds, since spin-up and spin-down orbitals 
are always equivalent due to magnetization density being 



identically zero (so that the exchange-correlation func- 
tional is spin non-polarized). Moreover, Eq. (|A8[) is iden- 
tical to Eq. (|A4[) . i.e., the densities of the m z = triplet 
and singlet states, when orbital relaxation effects are ig- 
nored. Since 7iNSp( r ) already provides a direct calcula- 
tion of an ensemble density, unlike the spin-purification 
formula which relies on a spurious mixed-state $mi,2, it 
can be expected that nNSp(r) could provide more accu- 
rate results. However, nNSp( r ) could in principle be a 
linear combination of other types of excited single deter- 
minant states. For instance, consider a doubly excited 
state <!>*, which is schematically represented in Fig. [TU 
The single-particle density tinsp ( r ) can be obtained from 



-\Lf 

-v., 



3> 



FIG. 14: Doubly excited-state that can participate in 
"•Nsp(r). 

the following ensemble density matrix 



5I*g)(*g| 



\ l*.X*.| 



(A9) 



More generally, jt-nsp (r) can be obtained from an ensem- 
ble of the form 



£> = ^£>l,0 + fA),0 



I|$ G )($ G | + i|$,)($, 



(A10) 

where q < 1. The reason why tinspO) results in better 
singlet energies could be related to the fact that it corre- 
sponds to an ensemble Eq. (|Al0j) . rather than a single de- 
terminant description of mixed spin states. Moreover, it 
is the only variational minimum that one can obtain using 
approximate ground-state exchange-correlation function- 
al which has this property. Due to the variational prin- 
ciple, one could expect that contributions from to be 
suppressed, since it is not the lowest energy configuration 
with vanishing spin polarization. However, it is not pos- 
sible to have control over the ensemble weights appearing 
in Eq. (|A10[) with approximate ground-state exchange- 
correlation functionals. Indeed, hints of this problem can 
be seen in Fig. [TT] where the singlet energy was found to 
increase with Ud almost linearly. One expects that larger 
Ud to destabilize the HOMO level, leading to a decreasing 
singlet energy. Due to the same reason, the ground-state 
energy should increase with Ud- Thus, an increasing sin- 
glet energy in Fig.[TT]is an indication that the NSP state 
contains ground-state contributions. The arbitrariness in 
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the ensembles that represent nNSp( r )> limits the under- tivatcd by its empirical success presented in section ITVl 
standing of the NSP A-SCF calculation. Therefore, the and previous studies in the literature [59M6l1 ] . 
use of riNSp( r ) is only partially justified, and mostly mo- 
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